%问题二求解：
data1 = readmatrix('附件','Range','A2:B1746');
ST=[9,10.5,12,13.5,15];
D=[306,337,0,31,61,92,122,153,184,214,245,275];
N=size(data1, 1);
lw=6;
l=lw;w=lw;
x0=0;y0=0;
xyz=[data1,zeros(N,1)+4];
e1=zeros(12,5);
% e1=0.9952;
e2=zeros(12,1);
e3=zeros(12,1);
e4=zeros(12,1);
e5=0.92;
for i=1:1
    for j=1:12
        [A,B]=SUN(10.5,D(j));
e1(j,i)=f1(xyz,l,w,N,A,B,x0,y0);%遮挡效率
e2(j,i)=f2(xyz,N,x0,y0,A,B);%余弦效率
e3(j,i)=f3(xyz,N,x0,y0);%大气透射率
e4(j,i)=f4(xyz,l,w,N,A,B,x0,y0);%截断效率
end
end
Se1=mean(e1);
Se2=mean(e2);
Se3=mean(e3);
Se4=mean(e4);
e=mean(Se2)*mean(Se3)*mean(Se4)*e1*e5;
W=e*1745*l*w;
data=data_func();
e11=zeros(12,5);
% e1=0.9952;
e22=zeros(12,1);
e33=zeros(12,1);
e44=zeros(12,1);
e55=0.92;
for i=1:1
    for j=1:12
        [A,B]=SUN(10.5,D(j));
e11(j,i)=f1(data,l,w,N,A,B,x0,y0);%遮挡效率
e22(j,i)=f2(data,N,x0,y0,A,B);%余弦效率
e33(j,i)=f3(data,N,x0,y0);%大气透射率
e44(j,i)=f4(data,l,w,N,A,B,x0,y0);%截断效率
end
end
ep2=mean(Se2)*mean(Se3)*mean(Se4)*e1*e5;

%fprintf("太阳光效率：%f\n",W);

% %问题一年平均功率
% e=readmatrix('效率','Range','G2:G13').*readmatrix('效率','Range','G15:G26').*readmatrix('效率','Range','G41:G52')*0.965015*0.92;
 ey=mean(ep2);
ST=[9,10.5,12,13.5,15];
D=[306,337,0,31,61,92,122,153,184,214,245,275];
DNI=zeros(12,5);
a=0.4237-0.00821*(6-3)^2;
b=0.5055+0.00595*(6.5-3)^2;
c=0.2711+0.01858*(2.5-3)^2;
for i=1:5
for j=1:12
[A,B]=SUN(ST(i),D(j));
DNI(j,i)=1.366*(a+b*exp(-c/sin(d(A))));
end
end
S=mean(DNI,2);
mean(DNI,'all');

aa=5.5.*5.5.*ones(12,5,size(data,1));
E=DNI.*aa.*e;%每块热功率，12*5*size(location,1)
Ef=sum(E,3);%镜场总瞬时热功率，12*5
year_Ef=sum(Ef,'all')/(12*5);%镜场年平均热功率,1*1
tmp=sum(aa,3);Ef_per_area=Ef/tmp(1);%单位面积输出热功率,12*5
%每块镜子年平均热功率，1745*1
mirror_year_ave=zeros(size(data,1),1);
% for i=1:size(data1,1)
% mirror_year_ave(i)=sum(E(:,:,i),'all')/(12*5);
% end
%单位面积输出月、年平均
month_Ef_per_area=zeros(12,1);
for i=1:12
tmp=Ef_per_area(i,:,:);
month_Ef_per_area(i)=sum(tmp,'all')/numel(tmp);
end
year_Ef_per_area=mean(month_Ef_per_area);


function [e4,SS]=f4(xyz,l,w,N,A,B,x0,y0)
%截断效率
z0=84;
a=[sind(B)*cosd(A),cosd(B)*cosd(A),sind(A)]; %入射光线
SS=zeros(N,1);
for i=1:N
sum1=0;sum2=0;
m=sqrt(xyz(i,1)^2+xyz(i,1)^2+xyz(i,3));
r=[-xyz(i,1)+x0,-xyz(i,2)+y0,-xyz(i,3)+z0]./sqrt((xyz(i,1)-x0)^2+(xyz(i,2)-y0)^2+(xyz(i,3)-z0)^2);%反射向量
n=real((r-a)/norm(r-a));%法向量
%镜面高度角、方位角
An=acosd(n(3));
Bn=atand(n(1)/n(2));
ss=[cosd(Bn),sind(Bn)*sind(An),-sind(Bn)*cosd(An);-sind(Bn),cosd(Bn)*sind(An),-cosd(Bn)*cosd(An);0,cosd(An),sind(An);]; %旋转矩阵
sss=[cosd(B),sind(B)*sind(A),sind(B)*cosd(A);-sind(B),cosd(B)*sind(A),cosd(B)*cosd(A);0,cosd(A),sind(A);]; %旋转矩阵
for j=1:10000 %模拟次数
AA=[sss*[rand(1)*4.64998*10^-3,rand(1)*2.16223*10^-5,1-rand(1)*2*10^-4]']';
R=2*dot(AA,n).*n-AA;
p=[ss*[-l*0.5*rand(1),-w*0.5*rand(1),0]']'+xyz(i,:);%随机一点
% 定义直线的参数方程
d=norm(cross([p(1)-x0,p(2)-y0,p(3)-z0],R))/sqrt(R(1)^2+R(2)^2+R(3)^2);
if d < 4
% 直线与圆柱曲面不相交';
sum1=sum1+1;
else
sum2=sum2+1;
end
end
SS(i)=sum1/10000;
end
e4=sum(SS)/N;
end


function [e3,SS]=f3(xyz,N,x0,y0)
%大气透射率
z0=84;
SS=zeros(length(xyz),1);
for i=1:N
dHR=sqrt((xyz(i,1)-x0)^2+(xyz(i,2)-y0)^2+(xyz(i,3)-z0)^2);
SS(i)=0.99321-0.0001176*dHR+1.97*10^-8*dHR^2;
end
e3=mean(SS);
end


function [e2,SS]=f2(xyz,N,x0,y0,A,B)
%余弦效率
z0=84;
SS=zeros(length(xyz),1);
for i=1:N
b=[xyz(i,1)-x0,xyz(i,2)-y0,xyz(i,3)-z0]; %反射光线
b=-1.*b;
a=[sind(B)*cosd(A),cosd(B)*cosd(A),sind(A)]; %入射光线
ta = acosd(dot(a,b)/(norm(a)*norm(b)));%入射光线与反射光线夹角
SS(i)=real(cosd(ta/2));
end
e2=mean(SS);
end

function e1=f1(xyz,l,w,N,A,B,x0,y0)
%遮挡效率
z0=84;
a=[sind(B)*cosd(A),cosd(B)*cosd(A),sind(A)]; %入射光
tt=zeros(length(xyz),5);
for i=1:size(xyz, 1)
    m=sqrt(xyz(i,1)^2+xyz(i,1)^2+xyz(i,3));
    r=[-xyz(i,1)+x0,-xyz(i,2)+y0,-xyz(i,3)+z0]/sqrt((xyz(i,1)-x0)^2+(xyz(i,2)-y0)^2+(xyz(i,3)-z0)^2);%反射向量
    n=real((r-a)/norm(r-a));%法向量
    %镜面高度角、方位角
    An=acos(n(3));
    Bn=atan(n(1)/n(2));
    ss=[cosd(Bn),sind(Bn)*sind(An),-sind(Bn)*cosd(An);-sind(Bn),cosd(Bn)*sind(An),-cosd(Bn)*cosd(An);0,cosd(An),sind(An);]; %翻转矩阵
    v1=[ss*[-0.5*l,-0.5*w,0]']'+xyz(i,:);%左下角
    v2=[ss*[0.5*l,-0.5*w,0]']'+xyz(i,:);
    v3=[ss*[-0.5*l,0.5*w,0]']'+xyz(i,:);
    tt(i,1)=v1(1);
    tt(i,2)=v1(2);
    tt(i,4)=real(abs(sum((v1-v2).*a)));
    tt(i,3)=real(abs(sum((v1-v3).*a)));
end
num = size(tt, 1);
% 计算重叠面积
for i = 1:num
    areas=0;
    for j = i+1:num
    area = rectint(tt(i,:), tt(j,:)); %两块之间的重叠
    areas = areas + area; %累加
    end
    tt(i,5)=areas;
end
SSS=(88/tand(A)-100)*7;
if SSS<0
    SSS=0;
end
e1=1-(sum(tt(:,5))+SSS)/N/l/w;
end
function f = canshu1(x)
%目标函数
N=2739;
H=4;
lw=6;
l=6;w=6;
x0=0;y0=-250;
ST=[9,10.5,12,13.5,15];
D=[306,337,0,31,61,92,122,153,184,214,245,275];
% for ii=1:N
xyz=[x(1:N)',x(N+1:2*N)',zeros(N,1)+x(3)];
% e1=0.9952;
e2=zeros(12,5);
% e3=zeros(12,5);
% e4=zeros(12,5);
% e5=0.92;
for i=1:5
    for j=1:12
        [A,B]=SUN(ST(i),D(j));
        e2(j,i)=f2(xyz,N,x0,y0,A,B);
    end
end
Se2=mean(e2,'all');
f=Se2;
end
function [A,B]=SUN(a,b)%B方位角，A高度角
        B1 = calculateSolarDeclination(b);%太阳赤纬角
        time = (pi/12.*(a-12));
        A1 = rad2deg(time);%太阳时角
        A = calculateHigh(A1,B1);%太阳高度角
        B = real(rad2deg((acos((sin(B1))-sin(A))*sin(deg2rad(39.4)))/(cos(A)*cos(deg2rad(39.4)))));

end
function delta = calculateSolarDeclination(element)
    epsilon = (23.45); 
    sinDelta = sin(2 .* pi.* element / 365) .* sin(2.*pi.*(epsilon)/360);
    delta = rad2deg(asin(sinDelta)) ;% 将结果转换为度
end
function x = d(aa)
    x = deg2rad(aa);
end
function High = calculateHigh(shijian,F)
linshi = deg2rad(39.4);
   bb = deg2rad(F);
   aa = deg2rad(shijian);
   High1 = cos(bb).*cos(linshi).*cos(aa)+sin(bb).*sin(linshi);
 High =rad2deg(asin(High1));
end

function data=data_func()% 这里求出的定日镜坐标是相对坐标
global W_a;
global L0;
R0 = 100;
delta_D=5+W_a+0.05;%给点容差%% !
coordinates = []; % 保存所有坐标
maxNumMirrors = 10000; % 估计的最大定日镜数
coordinates = NaN(maxNumMirrors, 2); % 预先分配空间
idx = 1; % 当前填充到的坐标索引
R=R0;
i=1;
while true
R = R + delta_D;
% 使用余弦定理计算 theta_D
 cos_theta_D = (2*R^2 - delta_D^2) / (2*R^2);
 theta_D = acos(cos_theta_D);
 ND = floor(2*pi/theta_D);
% % 重新计算 theta_D 的值
 theta_D = 2*pi/ND;
% % 计算 cos 和 sin 值
 cos_values = cos((0:ND-1) * theta_D);
 sin_values = sin((0:ND-1) * theta_D);
% % 计算所有定日镜的坐标
 x_values = R * cos_values;
 y_values = R * sin_values;
% % 更新坐标数组
 coordinates(idx:idx+ND-1, :) = [x_values', y_values'];
 idx = idx + ND;
%i = 1;
if R>(350+L0)
        break;
end 
distance_first_layer = sqrt((x_values(2)-x_values(1)^2+(y_values(2)-y_values(1))^2));
while true
prev=sqrt((coordinates(idx-1,1) - coordinates(idx-2,1))^2 + (coordinates(idx-1,2)-coordinates(idx-2,2))^2);
R = R + sqrt(delta_D^2 - (prev/2)^2);
if R > (350+L0)%%%%
break;
end
 theta_D0 = atan2(coordinates(idx-ND,2), coordinates(idx-ND,1)) + theta_D/2;
 cos_values = cos((0:ND-1) * theta_D + theta_D0); % 角平分线
 sin_values = sin((0:ND-1) * theta_D + theta_D0); % 角平分线
 x_values = R * cos_values;
 y_values = R * sin_values;
 coordinates(idx:idx+ND-1, :) = [x_values', y_values'];
 idx = idx + ND;
%distance_first_layer = sqrt((x_values(2)-x_values(1)^2+(y_values(2)-y_values(1))^2));
% 检查距离条件
distance = sqrt((x_values(2) - x_values(1))^2 + (y_values(2) - y_values(1))^2);
if distance > 1.33 * distance_first_layer
break;
end
layer = layer + 1;
end
i = i + 1;
end
% 删除未使用的空间
coordinates = coordinates(~isnan(coordinates(:,1)), :);
%筛选定日镜坐标
% 检查约束条件
N = size(coordinates,1);
valid_indices = [];
for i = 1:N
x1 = coordinates(i,1);
y1 = coordinates(i,2);
if x1^2 + (y1-L0)^2 <= 350^2
valid_indices = [valid_indices; i];
end
end
data=[coordinates(valid_indices,1), coordinates(valid_indices,2)];
error2 = check_distance(data);
if ~error2
    error('定日镜距离不符合');
end
data(:,3) = sqrt(data(:,1).^2 + data(:,2).^2);
data = sortrows(data, 3);
end
function result = check_distance(coordinates)% 检验点是否满足两个定日镜之间的距离要大于等于 delta_D
% 给定的参数
global W_a;
delta_D = 5 + W_a;
% 计算所有点之间的距离
distances = pdist(coordinates);
% 找到最小距离
min_distance = min(distances);
% 检查最小距离是否大于等于 delta_D
if min_distance >= delta_D
result = true;
else
    result = false;
end
end






